Having matched RNA and ATAC-seq profiles provides an opportunity to characterise the regulatory mechanisms that control somitogenesis.

## expression data
meta.rna <- read.table(paste0(dir, "RNA-seq/data/metadata_RNAseq.tsv"), 
                       stringsAsFactors = FALSE, header = TRUE)
meta.rna <- meta.rna[meta.rna$QC == 1 & meta.rna$wrongStage == 0,]
meta.rna$stage <- paste0("stage", meta.rna$stage)
meta.rna$somite <- factor(meta.rna$somite, levels=c("SIII","SII","SI"))

data.rna <- read.table(paste0(dir, "RNA-seq/data/geneCounts.NORM_batchCorrected_14PCs.tsv"),
                       check.names = FALSE)
data.rna <- data.rna[,c('gene', meta.rna$sample)]
universe <- data.rna$gene
names(universe) <- row.names(data.rna)

## accessibility data
meta.atac <- read.table(paste0(dir, "ATAC-seq/data/metadata_ATACseq_GQ.tsv"), 
                   stringsAsFactors = FALSE, header = TRUE)
meta.atac <- meta.atac[meta.atac$QCpass==1,]
meta.atac$stage <- paste0("stage", meta.atac$stage)
meta.atac$somite <- factor(meta.atac$somite, levels=c("SIII","SII","SI"))

data.atac <- readRDS(paste0(dir, "ATAC-seq/results/04_peakCounts_csawMerged.NORM.batchCorrected_18PCs.Rds"))

Peak to gene linking

One way to do this is to identify chromatin loci in the vicinity of DE genes that show coordinated activity patterns; these loci can be postulated as putative regulatory elements driving the changes in gene expression. We can exploit the variation across development, which usually involves large changes in activity, to find gene-peak pairs using correlation analysis. To determine whether a given correlation score is significantly larger than expected by chance, we compare correlation values of the same gene against peaks that are in other chromosomes, which do not contribute to cis regulatory mechanisms.

For this analysis, we use the method implemented in FigR (with minor modifications), which uses chromVAR to determine a set of background peaks matched for accessibility and GC content to model the null distribution of correlation scores for each gene-peak pair. Naturally, the analysis is restricted to the 43 samples where we generated good-quality libraries for both the RNA and ATAC datasets. We look for correlated peaks within 100kb of a DE gene.

There are 15734 gene-peak pairs that are significantly correlated. These involve nearly 6K DE genes, linked to ~13.5K peaks.

c(n_genes = length(unique(links$Gene)),
  n_peaks = length(unique(links$PeakRanges)))
## n_genes n_peaks 
##    5909   13657

Most gene-peak pairs show moderate correlation values (median=0.382547). However, a few pairs with very low correlation coefficients are deemed significant by the method. These correspond to cases where gene expression and accessibility are not associated (correlation values close to 0) but the set of background peaks used to compute significance are biased towards negative correlations. This leads to a significant p-value even if the correlation of the pair in question is not significant itself.

ggplot(links, aes(1, rObs)) +
  geom_violin() +
  geom_boxplot(width=0.1) +
  xlab("") +
  ylab("Spearman correlation") +
  th + theme(axis.text.x = element_blank(),
             axis.ticks.x = element_blank())

Thus, to consider a gene-peak pair as significantly correlated we require a minimum correlation coefficient of 0.3 on top of a significant p-value, taking the number of links to 12.8K.

links <- links[links$rObs>0.3,]
nrow(links)
## [1] 12803

Although most genes are linked to a single or a few peaks, a small proportion are connected to more than 5 peaks (349 genes). In contrast, the vast majority of peaks (88.91%) regulate a single gene.

par(mfrow=c(1,2))
plot(table(table(links$Gene)), bty="l",
     xlab="number of peaks per gene",
     ylab="number of genes")
plot(table(table(links$PeakRanges)), bty="l",
     xlab="number of genes per peak",
     ylab="number of peaks")

# summary(as.numeric(table(links$Gene)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 1.00    1.00    2.00    2.452   3.000  30.000  
# summary(as.numeric(table(links$PeakRanges)))
# Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 1.000   1.000   1.000   1.147   1.000   8.000 

The set of defined pairs involve around half of all the DE genes, for both those that are differential across somite trios, and across stages.

## check how many of the correlated regions are differential
dars.somite <- read.table(paste0(dir, "RNA+ATAC/results/01_DAregions_summary_somiteTrios_withClusters.tsv"))
dars.stage <- read.table(paste0(dir, "RNA+ATAC/results/02_DAregions_summary_stage_fate.tsv"))
dars <- union(row.names(dars.somite), row.names(dars.stage))

links$DE <- ifelse(links$Gene %in% degs.somite$gene & links$Gene %in% degs.stage$gene, "both",
                   ifelse(links$Gene %in% degs.somite$gene, "somite",
                          ifelse(links$Gene %in% degs.stage$gene, "stage", "nonDE")))
links$DE_stage <- degs.stage[match(links$Gene, degs.stage$gene),]$fate 
# the same proportion of DE genes in each fate are linked

links$DA <- ifelse(links$PeakRanges %in% row.names(dars.somite) & links$PeakRanges %in% row.names(dars.stage), "both",
                   ifelse(links$PeakRanges %in% row.names(dars.somite), "somite",
                          ifelse(links$PeakRanges %in% row.names(dars.stage), "stage", "nonDA")))
links$DA_stage <- dars.stage[match(links$PeakRanges, row.names(dars.stage)),]$fate 

df <- data.frame(n_DE_linked = c(length(unique(links[links$DE %in% c("both", "somite"),]$Gene)),
                                 length(unique(links[links$DE %in% c("both", "stage"),]$Gene))),
                 DE_pct = c(round(length(unique(links[links$DE %in% c("both", "somite"),]$Gene))/nrow(degs.somite)*100,2),
                            round(length(unique(links[links$DE %in% c("both", "stage"),]$Gene))/nrow(degs.stage)*100,2)),
                 n_DA_linked = c(length(unique(links[links$DA %in% c("both", "somite"),]$PeakRanges)),
                                 length(unique(links[links$DA %in% c("both", "stage"),]$PeakRanges))),
                 DA_pct = c(round(length(unique(links[links$DA %in% c("both", "somite"),]$PeakRanges))/nrow(dars.somite)*100,2),
                            round(length(unique(links[links$DA %in% c("both", "stage"),]$PeakRanges))/nrow(dars.stage)*100,2)))
row.names(df) <- c("somite", "stage")
df

Although some peaks are nearby their putative target, the majority are dozens of kilobases away.

## TSS coordinates
TSSg <- FigR::mm10TSSRanges

# distance between gene and peak
tmp <- unlist(lapply(strsplit(links$PeakRanges, ":"), '[[', 2))
links.gr <- GenomicInteractions(GRanges(seqnames = unlist(lapply(strsplit(links$PeakRanges, ":"), '[[', 1)),
                                        IRanges(start = as.numeric(unlist(lapply(strsplit(tmp, "-"), '[[', 1))),
                                                end = as.numeric(unlist(lapply(strsplit(tmp, "-"), '[[', 2))))),
                                TSSg[match(links$Gene, TSSg$gene_name),])
distance <- distance(anchorOne(links.gr), anchorTwo(links.gr))
links$distance <- distance

bins <- data.frame(distance=c('0-200','200-2e3','2e3-5e3','5e3-10e3','10e3-50e3','50e3-100e3'),
                   n=c(sum(distance <= 200),
                       sum(distance > 200 & distance <= 2e3),
                       sum(distance > 2e3 & distance <= 5e3),
                       sum(distance > 5e3 & distance <= 10e3),
                       sum(distance > 10e3 & distance <= 50e3),
                       sum(distance > 50e3 & distance <= 100e3)))
bins$n_pct <- round(bins$n/length(distance)*100,2)
bins

The median distance between linked peaks and genes is of 37.845kb.

Supporting the regulatory potential of the linked peaks, a larger proportion overlap ENCODE cCREs characterised as enhancers (H3K27ac-high, H3K4me3-low) compared to all peaks detected in somites, and all peaks within 100kb of a DE gene. Furthermore, the overlap greatly increases if we restrict to links with the strongest correlations.

encode <- read.table(paste0(dir, "ATAC-seq/results/05_peaks_overlapExternalData.tsv"), sep="\t")
row.names(encode) <- paste0("chr", row.names(encode))

encode$encode <- ifelse(encode$pELS | encode$dELS, "enhancer",
                        ifelse(encode$PLS, "promoter",
                               ifelse(encode$CTCF_bound, "CTCFonly",
                                      ifelse(encode$H3K4me3, "H3K4me3only", "not_in_encode"))))
links$encode <- encode[match(links$PeakRanges, row.names(encode)),]$encode
links$fantom <- encode[match(links$PeakRanges, row.names(encode)),]$FANTOM
links$vista <- encode[match(links$PeakRanges, row.names(encode)),]$VISTA

links$encode <- as.factor(links$encode)

# round(rbind(all_peaks = prop.table(table(encode$encode)[1:4])*100,
#             peaks_proximal_to_DEGs = prop.table(table(encode[unique(cisCor$PeakRanges),]$encode)[1:4])*100,
#             peaks_linked_to_DEGs = prop.table(table(links$encode)[1:4])*100,
#             linked_0.5 = prop.table(table(links[links$rObs>0.5,]$encode)[1:4])*100,
#             linked_0.65 = prop.table(table(links[links$rObs>0.65,]$encode)[1:4])*100), 2)[,c(4,2,1,3)]

df <- data.frame(type = rep(c("enhancer", "not_in_encode"), each=7),
                 prop = c(sum(encode$encode=="enhancer")/nrow(encode)*100,
                          sum(encode[unique(cisCor$PeakRanges),]$encode=="enhancer")/nrow(encode[unique(cisCor$PeakRanges),])*100,
                          sum(links[links$rObs>0.3,]$encode=="enhancer")/nrow(links[links$rObs>0.3,])*100,
                          sum(links[links$rObs>0.4,]$encode=="enhancer")/nrow(links[links$rObs>0.4,])*100,
                          sum(links[links$rObs>0.5,]$encode=="enhancer")/nrow(links[links$rObs>0.5,])*100,
                          sum(links[links$rObs>0.6,]$encode=="enhancer")/nrow(links[links$rObs>0.6,])*100,
                          sum(links[links$rObs>0.7,]$encode=="enhancer")/nrow(links[links$rObs>0.7,])*100,
                          sum(encode$encode=="not_in_encode")/nrow(encode)*100,
                          sum(encode[unique(cisCor$PeakRanges),]$encode=="not_in_encode")/nrow(encode[unique(cisCor$PeakRanges),])*100,
                          sum(links[links$rObs>0.3,]$encode=="not_in_encode")/nrow(links[links$rObs>0.3,])*100,
                          sum(links[links$rObs>0.4,]$encode=="not_in_encode")/nrow(links[links$rObs>0.4,])*100,
                          sum(links[links$rObs>0.5,]$encode=="not_in_encode")/nrow(links[links$rObs>0.5,])*100,
                          sum(links[links$rObs>0.6,]$encode=="not_in_encode")/nrow(links[links$rObs>0.6,])*100,
                          sum(links[links$rObs>0.7,]$encode=="not_in_encode")/nrow(links[links$rObs>0.7,])*100))

ggplot(df, aes(rep(1:7,2), prop, colour=type)) +
  geom_point(size=2) +
  geom_line() +
  xlab("") +
  ylab("% of peaks") +
  ylim(c(0,100)) +
  scale_x_continuous(breaks=1:7,
                     labels=c("all peaks", "peaks near DEGs", paste0("linked peaks > 0.", 3:7))) +
  th + theme(axis.text.x = element_text(angle=45, hjust=1),
             panel.grid.major.y = element_line(size=0.5, colour="grey90"),
             panel.grid.minor.y = element_line(size=0.25, colour="grey95"))

Similarly, while only ~10% of all peaks overlap with FANTOM5 enhancers, the number increases to 14% for linked peaks and this goes up as the correlation threshold becomes more stringent.

df <- data.frame(type = "in_FANTOM",
                 prop = c(sum(encode$FANTOM)/nrow(encode)*100,
                          sum(encode[unique(cisCor$PeakRanges),]$FANTOM)/nrow(encode[unique(cisCor$PeakRanges),])*100,
                          sum(links[links$rObs>0.3,]$fantom)/nrow(links[links$rObs>0.3,])*100,
                          sum(links[links$rObs>0.4,]$fantom)/nrow(links[links$rObs>0.4,])*100,
                          sum(links[links$rObs>0.5,]$fantom)/nrow(links[links$rObs>0.5,])*100,
                          sum(links[links$rObs>0.6,]$fantom)/nrow(links[links$rObs>0.6,])*100,
                          sum(links[links$rObs>0.7,]$fantom)/nrow(links[links$rObs>0.7,])*100))

ggplot(df, aes(1:7, prop, colour=type)) +
  geom_point(size=2) +
  geom_line() +
  xlab("") +
  ylab("% of peaks") +
  ylim(c(0,100)) +
  scale_x_continuous(breaks=1:7,
                     labels=c("all peaks", "peaks near DEGs", paste0("linked peaks > 0.", 3:7))) +
  th + theme(axis.text.x = element_text(angle=45, hjust=1),
             panel.grid.major.y = element_line(size=0.5, colour="grey90"),
             panel.grid.minor.y = element_line(size=0.25, colour="grey95"))

Highly-regulated genes

As mentioned above, although most genes are linked only to a single or few peaks, a few hundred have many more connections. The authors of FigR coined the term domains of regulatory chromatin (DORCs) to refer to these loci of high interactivity, and showed that they “are enriched for lineage-determining genes and overlap with known super-enhancers”. In our case, many of the genes linked to dozens of peaks correspond to late Hox genes, which is expected since it has been shown that Hox expression regulation relies on coordinated chromatin remodeling. Many other genes are identified as well.

# Determine DORC genes
dorcGenes <- dorcJPlot(dorcTab = links,
                       cutoff = 6,
                       returnGeneList = TRUE)

This restricted set of genes with a large number of linked peaks are still enriched for key processes relevant for somitogenesis, such as patterning and the development and morphogenesis of the different somite derivative lineages. These data suggest that the expression of genes that are crucial in the development and differentiation of somites are under intricate regulatory control.

go <- topGOtable(DEgenes = dorcGenes,
                 BGgenes = universe,
                 topGO_method2 = "elim",
                 ontology = "BP",
                 geneID = "symbol",
                 fullNamesInRows = FALSE,
                 addGeneToTerms = TRUE,
                 mapping = "org.Mm.eg.db",
                 topTablerows = 100)
go[,c(2,4,7,9)]

Interestingly, while genes most active at different fates are all as likely to be associated to at least one peak, there is a strong depletion of genes expressed in thoracic somites among DORC genes.

## check proportion per stage
tmp <- unique(links[,c('Gene','DE_stage')])
df <- rbind(all_DEGs = prop.table(table(degs.stage$fate))*100,
            all_linked_DEGs = prop.table(table(tmp$DE_stage))*100,
            DORC_DEGs  = prop.table(table(tmp[tmp$Gene %in% dorcGenes,]$DE_stage))*100)
df[,c(1,4,2,3)]
##                 cervical  thoracic   lumbar   sacral
## all_DEGs        34.49304 17.902584 17.58449 30.01988
## all_linked_DEGs 32.70559 15.861059 17.99540 33.43796
## DORC_DEGs       38.05310  8.259587 15.63422 38.05310

Regulatory control of DORC genes

Next, we can ask whether the peaks linked to the highly-regulated genes are enriched for particular TF binding sites, compared to a set of (matched) background peaks. These TFs have the potential of regulating the expression of the gene in question. This hypothesis becomes stronger if the expression of the TF itself is correlated with the accessibility values of the peaks. FigR performs these two calculations, and combines the results into a regulation score. Positive/negative scores correspond to activating/repressing interactions.

We observe a set of DORCs with very strong TFBS enrichment scores, all corresponding to enrichment for Nr6a1 binding sites; DORCs regulated by Nr6a1 show negative scores, indicating TF expression is negatively correlated with chromatin accessibility. Other transcription factors result in more modest motif enrichment scores, with both repressing and enhancing activities.

ggplot(fig.d, aes(Corr.log10P, Enrichment.log10P, color=Score, shape=as.factor(Motif=="Nr6a1"))) +
  geom_point_rast(size=0.75) + 
  scale_color_gradientn(colours = jdb_palette("solar_extra"), 
                        limits=c(-2,2),
                        oob = scales::squish,
                        breaks=scales::breaks_pretty(n=3)) +
  geom_vline(xintercept = 0, lty=2, colour="grey") +
  xlab(expression('(signed) -log'[10]*' p-value TF expression - accessibility correlation')) +
  ylab(expression('-log'[10]*' p-value TF motif enrichment')) +
  labs(shape = "Nr6a1 enrichment") +
  th

Across all DORC genes, some TFs consistently show large regulation scores, as shown below, with Nr6a1 again standing out. Several of the TFs highlighted have been already picked up in previous analyses, given the enrichment of their binding sites in DA peaks.

rankDrivers(fig.d, rankBy = "meanScore", interactive = FALSE)

Generally, TFs are associated with a handful to a dozen different targets, but some factors show more widespread activity (again, Nr6a1, Ebf1 and Sox6, all acting as repressors, and several Hox genes).

rankDrivers(fig.d, score.cut = 1, rankBy = "nTargets", interactive = TRUE)
## Ranking TFs by total number of associated DORCs ..
## Using absolute score cut-off of: 1 ..
## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

We can then group DORCs (rows) that are regulated in a concerted manner by groups of TFs (columns).

score.cut <- 1.25
set.seed(2948)
regulons <- plotfigRHeatmap(figR.d = fig.d,
                            score.cut = score.cut,
                            row_names_gp = gpar(fontsize=8, fontface="italic"),
                            column_names_gp = gpar(fontsize=10),
                            show_row_dend = TRUE, k=4)
regulons <- draw(regulons)

## get cluster membership
DORCsToKeep <- fig.d %>% dplyr::filter(abs(Score) >= score.cut) %>% pull(DORC) %>% unique()
TFsToKeep <- fig.d %>% dplyr::filter(abs(Score) >= score.cut) %>% pull(Motif) %>% unique()
net.d <- fig.d %>% dplyr::filter(DORC %in% DORCsToKeep & Motif %in% TFsToKeep) %>%
    reshape2::dcast(DORC ~ Motif) %>%
    tibble::column_to_rownames("DORC") %>% as.matrix()

clusters <- as.data.frame(unlist(row_order(regulons)))
clusters$cluster <- paste0("cluster", substr(row.names(clusters), 1, 1))
clusters$gene <- row.names(net.d[clusters[,1],])
clusters$cluster <- factor(clusters$cluster, levels=paste0("cluster", names(row_order(regulons))))

And visualise the expression of the DORC genes, to link each module to its activity pattern across stages.

data <- data.rna[match(clusters$gene, row.names(data.rna)),]
data <- t(apply(data, 1, function(x) (x-mean(x))/sd(x) ))

m <- meta.rna[match(colnames(data), meta.rna$sample),]
m$stage <- factor(m$stage, levels=paste0("stage", c(8,18,21,25,27,35)))
m$somite <- factor(m$somite, levels=c("SIII","SII","SI"))
order <- order(m$stage, m$somite)

## heatmap annotation
ha  <- HeatmapAnnotation(df = data.frame(stage = m[order,]$stage, 
                                         somite = m[order,]$somite), 
                         col = list(stage = cols.stage, somite = cols.somite))
Heatmap(data[,order], 
        cluster_columns = FALSE, 
        cluster_rows = FALSE,
        col=colorRamp2(breaks = c(-4,-2,0,2,3), 
                       colors = c("steelblue4","steelblue","white", 
                                  "indianred3","indianred4")),
        row_names_gp = gpar(fontsize=8, fontface="italic"),
        column_names_gp = gpar(fontsize=10),
        name = "z-score", 
        show_row_names = TRUE, 
        show_column_names = FALSE, 
        top_annotation = ha,
        row_split = clusters$cluster,
        cluster_row_slices = FALSE)

Two main modules are evident, with genes expressed either early or late in development. The role of Nr6a1 is prominent, and Hox genes also show regulatory effects consistent with developmental progression.

saveRDS(cisCor, paste0(dir, "RNA+ATAC/results/03_gene_peak_correlations.Rds"))
write.table(links, paste0(dir, "RNA+ATAC/results/03_gene_peak_links.tsv"), quote = FALSE, sep = "\t", row.names = FALSE)
saveRDS(fig.d, paste0(dir, "RNA+ATAC/results/03_FigR_network.Rds"))
saveRDS(net.d, paste0(dir, "RNA+ATAC/results/03_regulationScores_thr1.25.Rds"))
write.table(clusters, paste0(dir, "RNA+ATAC/results/03_regulons_clusters.tsv"), quote = FALSE, sep="\t", row.names = FALSE)
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] BSgenome.Mmusculus.UCSC.mm10_1.4.0 BSgenome_1.60.0                   
##  [3] rtracklayer_1.52.1                 Biostrings_2.60.2                 
##  [5] XVector_0.32.0                     org.Mm.eg.db_3.13.0               
##  [7] FigR_0.1.0                         motifmatchr_1.14.0                
##  [9] chromVAR_1.14.0                    cowplot_1.1.1                     
## [11] Matrix_1.4-1                       pbmcapply_1.5.0                   
## [13] doParallel_1.0.17                  iterators_1.0.14                  
## [15] foreach_1.5.2                      dynamicTreeCut_1.63-1             
## [17] pcaExplorer_2.18.0                 topGO_2.44.0                      
## [19] SparseM_1.81                       GO.db_3.13.0                      
## [21] AnnotationDbi_1.54.1               graph_1.70.0                      
## [23] networkD3_0.4                      circlize_0.4.14                   
## [25] ComplexHeatmap_2.8.0               RColorBrewer_1.1-3                
## [27] BuenColors_0.5.6                   MASS_7.3-56                       
## [29] GGally_2.1.2                       ggrastr_1.0.1                     
## [31] ggpubr_0.4.0                       ggrepel_0.9.1                     
## [33] ggplot2_3.3.5                      dplyr_1.0.8                       
## [35] GenomicInteractions_1.26.0         InteractionSet_1.20.0             
## [37] SummarizedExperiment_1.22.0        Biobase_2.52.0                    
## [39] GenomicRanges_1.44.0               GenomeInfoDb_1.28.4               
## [41] IRanges_2.26.0                     S4Vectors_0.30.2                  
## [43] BiocGenerics_0.38.0                MatrixGenerics_1.4.3              
## [45] matrixStats_0.62.0                
## 
## loaded via a namespace (and not attached):
##   [1] rappdirs_0.3.3              AnnotationForge_1.34.1     
##   [3] R.methodsS3_1.8.1           nabor_0.5.0                
##   [5] pkgmaker_0.32.2             tidyr_1.2.0                
##   [7] bit64_4.0.5                 knitr_1.38                 
##   [9] R.utils_2.11.0              DelayedArray_0.18.0        
##  [11] data.table_1.14.2           rpart_4.1.16               
##  [13] KEGGREST_1.32.0             TFBSTools_1.30.0           
##  [15] RCurl_1.98-1.6              AnnotationFilter_1.16.0    
##  [17] generics_0.1.2              snow_0.4-4                 
##  [19] GenomicFeatures_1.44.2      RSQLite_2.2.12             
##  [21] tzdb_0.3.0                  bit_4.0.4                  
##  [23] webshot_0.5.3               xml2_1.3.3                 
##  [25] httpuv_1.6.5                assertthat_0.2.1           
##  [27] DirichletMultinomial_1.34.0 viridis_0.6.2              
##  [29] xfun_0.30                   hms_1.1.1                  
##  [31] jquerylib_0.1.4             evaluate_0.15              
##  [33] promises_1.2.0.1            TSP_1.2-0                  
##  [35] fansi_1.0.3                 restfulr_0.0.13            
##  [37] progress_1.2.2              caTools_1.18.2             
##  [39] dendextend_1.15.2           dbplyr_2.1.1               
##  [41] Rgraphviz_2.36.0            igraph_1.3.0               
##  [43] DBI_1.1.2                   geneplotter_1.70.0         
##  [45] htmlwidgets_1.5.4           reshape_0.8.9              
##  [47] Rmpfr_0.8-7                 purrr_0.3.4                
##  [49] ellipsis_0.3.2              crosstalk_1.2.0            
##  [51] backports_1.4.1             annotate_1.70.0            
##  [53] gridBase_0.4-7              biomaRt_2.48.3             
##  [55] vctrs_0.4.1                 ensembldb_2.16.4           
##  [57] Cairo_1.5-15                abind_1.4-5                
##  [59] cachem_1.0.6                withr_2.5.0                
##  [61] Gviz_1.36.2                 checkmate_2.0.0            
##  [63] GenomicAlignments_1.28.0    prettyunits_1.1.1          
##  [65] cluster_2.1.3               seqLogo_1.58.0             
##  [67] lazyeval_0.2.2              crayon_1.5.1               
##  [69] genefilter_1.74.1           labeling_0.4.2             
##  [71] pkgconfig_2.0.3             vipor_0.4.5                
##  [73] ProtGenerics_1.24.0         seriation_1.3.5            
##  [75] nnet_7.3-17                 rlang_1.0.2                
##  [77] lifecycle_1.0.1             miniUI_0.1.1.1             
##  [79] registry_0.5-1              filelock_1.0.2             
##  [81] BiocFileCache_2.0.0         doSNOW_1.0.20              
##  [83] GOstats_2.58.0              dichromat_2.0-0            
##  [85] rngtools_1.5.2              carData_3.0-5              
##  [87] base64enc_0.1-3             beeswarm_0.4.0             
##  [89] GlobalOptions_0.1.2         pheatmap_1.0.12            
##  [91] png_0.1-7                   viridisLite_0.4.0          
##  [93] rjson_0.2.21                bitops_1.0-7               
##  [95] shinydashboard_0.7.2        R.oo_1.24.0                
##  [97] blob_1.2.3                  shape_1.4.6                
##  [99] stringr_1.4.0               readr_2.1.2                
## [101] jpeg_0.1-9                  rstatix_0.7.0              
## [103] shinyAce_0.4.1              ggsignif_0.6.3             
## [105] CNEr_1.28.0                 scales_1.2.0               
## [107] memoise_2.0.1               GSEABase_1.54.0            
## [109] magrittr_2.0.3              plyr_1.8.7                 
## [111] zlibbioc_1.38.0             threejs_0.3.3              
## [113] compiler_4.1.0              BiocIO_1.2.0               
## [115] clue_0.3-60                 DESeq2_1.32.0              
## [117] Rsamtools_2.8.0             cli_3.2.0                  
## [119] Category_2.58.0             htmlTable_2.4.0            
## [121] Formula_1.2-4               tidyselect_1.1.2           
## [123] stringi_1.7.6               shinyBS_0.61.1             
## [125] highr_0.9                   yaml_2.3.5                 
## [127] locfit_1.5-9.5              latticeExtra_0.6-29        
## [129] sass_0.4.1                  VariantAnnotation_1.38.0   
## [131] tools_4.1.0                 rstudioapi_0.13            
## [133] TFMPvalue_0.0.8             foreign_0.8-82             
## [135] gridExtra_2.3               farver_2.1.0               
## [137] digest_0.6.29               FNN_1.1.3                  
## [139] pracma_2.3.8                shiny_1.7.1                
## [141] Rcpp_1.0.8.3                car_3.0-12                 
## [143] broom_0.8.0                 later_1.3.0                
## [145] httr_1.4.2                  biovizBase_1.40.0          
## [147] colorspace_2.0-3            XML_3.99-0.9               
## [149] splines_4.1.0               RBGL_1.68.0                
## [151] plotly_4.10.0               gmp_0.6-5                  
## [153] xtable_1.8-4                poweRlaw_0.70.6            
## [155] jsonlite_1.8.0              heatmaply_1.3.0            
## [157] R6_2.5.1                    Hmisc_4.7-0                
## [159] pillar_1.7.0                htmltools_0.5.2            
## [161] mime_0.12                   NMF_0.24.0                 
## [163] glue_1.6.2                  fastmap_1.1.0              
## [165] DT_0.22                     BiocParallel_1.26.2        
## [167] codetools_0.2-18            utf8_1.2.2                 
## [169] lattice_0.20-45             bslib_0.3.1                
## [171] tibble_3.1.6                curl_4.3.2                 
## [173] ggbeeswarm_0.6.0            gtools_3.9.2               
## [175] magick_2.7.3                survival_3.3-1             
## [177] limma_3.48.3                rmarkdown_2.13             
## [179] munsell_0.5.0               GetoptLong_1.0.5           
## [181] GenomeInfoDbData_1.2.6      reshape2_1.4.4             
## [183] gtable_0.3.0